home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / dtgsen.z / dtgsen
Encoding:
Text File  |  2002-10-03  |  17.7 KB  |  463 lines
  1.  
  2.  
  3.  
  4. DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      DTGSEN - reorder the generalized real Schur decomposition of a real
  10.      matrix pair (A, B) (in terms of an orthonormal equivalence trans-
  11.      formation Q' * (A, B) * Z), so that a selected cluster of eigenvalues
  12.      appears in the leading diagonal blocks of the upper quasi-triangular
  13.      matrix A and the upper triangular B
  14.  
  15. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  16.      SUBROUTINE DTGSEN( IJOB, WANTQ, WANTZ, SELECT, N, A, LDA, B, LDB, ALPHAR,
  17.                         ALPHAI, BETA, Q, LDQ, Z, LDZ, M, PL, PR, DIF, WORK,
  18.                         LWORK, IWORK, LIWORK, INFO )
  19.  
  20.          LOGICAL        WANTQ, WANTZ
  21.  
  22.          INTEGER        IJOB, INFO, LDA, LDB, LDQ, LDZ, LIWORK, LWORK, M, N
  23.  
  24.          DOUBLE         PRECISION PL, PR
  25.  
  26.          LOGICAL        SELECT( * )
  27.  
  28.          INTEGER        IWORK( * )
  29.  
  30.          DOUBLE         PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), B(
  31.                         LDB, * ), BETA( * ), DIF( * ), Q( LDQ, * ), WORK( * ),
  32.                         Z( LDZ, * )
  33.  
  34. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  35.      These routines are part of the SCSL Scientific Library and can be loaded
  36.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  37.      directs the linker to use the multi-processor version of the library.
  38.  
  39.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  40.      4 bytes (32 bits). Another version of SCSL is available in which integers
  41.      are 8 bytes (64 bits).  This version allows the user access to larger
  42.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  43.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  44.      only one of the two versions; 4-byte integer and 8-byte integer library
  45.      calls cannot be mixed.
  46.  
  47. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  48.      DTGSEN reorders the generalized real Schur decomposition of a real matrix
  49.      pair (A, B) (in terms of an orthonormal equivalence trans- formation Q' *
  50.      (A, B) * Z), so that a selected cluster of eigenvalues appears in the
  51.      leading diagonal blocks of the upper quasi-triangular matrix A and the
  52.      upper triangular B. The leading columns of Q and Z form orthonormal bases
  53.      of the corresponding left and right eigen- spaces (deflating subspaces).
  54.      (A, B) must be in generalized real Schur canonical form (as returned by
  55.      DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 diagonal
  56.      blocks. B is upper triangular.
  57.  
  58.      DTGSEN also computes the generalized eigenvalues
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))
  71.  
  72.  
  73.  
  74.                  w(j) = (ALPHAR(j) + i*ALPHAI(j))/BETA(j)
  75.  
  76.      of the reordered matrix pair (A, B).
  77.  
  78.      Optionally, DTGSEN computes the estimates of reciprocal condition numbers
  79.      for eigenvalues and eigenspaces. These are Difu[(A11,B11), (A22,B22)] and
  80.      Difl[(A11,B11), (A22,B22)], i.e. the separation(s) between the matrix
  81.      pairs (A11, B11) and (A22,B22) that correspond to the selected cluster
  82.      and the eigenvalues outside the cluster, resp., and norms of
  83.      "projections" onto left and right eigenspaces w.r.t.  the selected
  84.      cluster in the (1,1)-block.
  85.  
  86.  
  87. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  88.      IJOB    (input) INTEGER
  89.              Specifies whether condition numbers are required for the cluster
  90.              of eigenvalues (PL and PR) or the deflating subspaces (Difu and
  91.              Difl):
  92.              =0: Only reorder w.r.t. SELECT. No extras.
  93.              =1: Reciprocal of norms of "projections" onto left and right
  94.              eigenspaces w.r.t. the selected cluster (PL and PR).  =2: Upper
  95.              bounds on Difu and Difl. F-norm-based estimate
  96.              (DIF(1:2)).
  97.              =3: Estimate of Difu and Difl. 1-norm-based estimate
  98.              (DIF(1:2)).  About 5 times as expensive as IJOB = 2.  =4: Compute
  99.              PL, PR and DIF (i.e. 0, 1 and 2 above): Economic version to get
  100.              it all.  =5: Compute PL, PR and DIF (i.e. 0, 1 and 3 above)
  101.  
  102.      WANTQ   (input) LOGICAL
  103.  
  104.      WANTZ   (input) LOGICAL
  105.  
  106.      SELECT  (input) LOGICAL array, dimension (N)
  107.              SELECT specifies the eigenvalues in the selected cluster.  To
  108.              select a real eigenvalue w(j), SELECT(j) must be set to w(j) and
  109.              w(j+1), corresponding to a 2-by-2 diagonal block, either
  110.              SELECT(j) or SELECT(j+1) or both must be set to either both
  111.              included in the cluster or both excluded.
  112.  
  113.      N       (input) INTEGER
  114.              The order of the matrices A and B. N >= 0.
  115.  
  116.      A       (input/output) DOUBLE PRECISION array, dimension(LDA,N)
  117.              On entry, the upper quasi-triangular matrix A, with (A, B) in
  118.              generalized real Schur canonical form.  On exit, A is overwritten
  119.              by the reordered matrix A.
  120.  
  121.      LDA     (input) INTEGER
  122.              The leading dimension of the array A. LDA >= max(1,N).
  123.  
  124.  
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))
  137.  
  138.  
  139.  
  140.      B       (input/output) DOUBLE PRECISION array, dimension(LDB,N)
  141.              On entry, the upper triangular matrix B, with (A, B) in
  142.              generalized real Schur canonical form.  On exit, B is overwritten
  143.              by the reordered matrix B.
  144.  
  145.      LDB     (input) INTEGER
  146.              The leading dimension of the array B. LDB >= max(1,N).
  147.  
  148.      ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
  149.              ALPHAI  (output) DOUBLE PRECISION array, dimension (N) BETA
  150.              (output) DOUBLE PRECISION array, dimension (N) On exit,
  151.              (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will be the
  152.              generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i and
  153.              BETA(j),j=1,...,N  are the diagonals of the complex Schur form
  154.              (S,T) that would result if the 2-by-2 diagonal blocks of the real
  155.              generalized Schur form of (A,B) were further reduced to
  156.              triangular form using complex unitary transformations.  If
  157.              ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive,
  158.              then the j-th and (j+1)-st eigenvalues are a complex conjugate
  159.              pair, with ALPHAI(j+1) negative.
  160.  
  161.      Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
  162.              On entry, if WANTQ = .TRUE., Q is an N-by-N matrix.  On exit, Q
  163.              has been postmultiplied by the left orthogonal transformation
  164.              matrix which reorder (A, B); The leading M columns of Q form
  165.              orthonormal bases for the specified pair of left eigenspaces
  166.              (deflating subspaces).  If WANTQ = .FALSE., Q is not referenced.
  167.  
  168.      LDQ     (input) INTEGER
  169.              The leading dimension of the array Q.  LDQ >= 1; and if WANTQ =
  170.              .TRUE., LDQ >= N.
  171.  
  172.      Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
  173.              On entry, if WANTZ = .TRUE., Z is an N-by-N matrix.  On exit, Z
  174.              has been postmultiplied by the left orthogonal transformation
  175.              matrix which reorder (A, B); The leading M columns of Z form
  176.              orthonormal bases for the specified pair of left eigenspaces
  177.              (deflating subspaces).  If WANTZ = .FALSE., Z is not referenced.
  178.  
  179.      LDZ     (input) INTEGER
  180.              The leading dimension of the array Z. LDZ >= 1; If WANTZ =
  181.              .TRUE., LDZ >= N.
  182.  
  183.      M       (output) INTEGER
  184.              The dimension of the specified pair of left and right eigen-
  185.              spaces (deflating subspaces). 0 <= M <= N.
  186.  
  187.              PL, PR  (output) DOUBLE PRECISION If IJOB = 1, 4 or 5, PL, PR are
  188.              lower bounds on the reciprocal of the norm of "projections" onto
  189.              left and right eigenspaces with respect to the selected cluster.
  190.              0 < PL, PR <= 1.  If M = 0 or M = N, PL = PR  = 1.  If IJOB = 0,
  191.              2 or 3, PL and PR are not referenced.
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))
  203.  
  204.  
  205.  
  206.      DIF     (output) DOUBLE PRECISION array, dimension (2).
  207.              If IJOB >= 2, DIF(1:2) store the estimates of Difu and Difl.
  208.              If IJOB = 2 or 4, DIF(1:2) are F-norm-based upper bounds on
  209.              Difu and Difl. If IJOB = 3 or 5, DIF(1:2) are 1-norm-based
  210.              estimates of Difu and Difl.  If M = 0 or N, DIF(1:2) = F-norm([A,
  211.              B]).  If IJOB = 0 or 1, DIF is not referenced.
  212.  
  213.      WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  214.              IF IJOB = 0, WORK is not referenced.  Otherwise, on exit, if INFO
  215.              = 0, WORK(1) returns the optimal LWORK.
  216.  
  217.      LWORK   (input) INTEGER
  218.              The dimension of the array WORK. LWORK >=  4*N+16.  If IJOB = 1,
  219.              2 or 4, LWORK >= MAX(4*N+16, 2*M*(N-M)).  If IJOB = 3 or 5, LWORK
  220.              >= MAX(4*N+16, 4*M*(N-M)).
  221.  
  222.              If LWORK = -1, then a workspace query is assumed; the routine
  223.              only calculates the optimal size of the WORK array, returns this
  224.              value as the first entry of the WORK array, and no error message
  225.              related to LWORK is issued by XERBLA.
  226.  
  227.      IWORK   (workspace/output) INTEGER array, dimension (LIWORK)
  228.              IF IJOB = 0, IWORK is not referenced.  Otherwise, on exit, if
  229.              INFO = 0, IWORK(1) returns the optimal LIWORK.
  230.  
  231.      LIWORK  (input) INTEGER
  232.              The dimension of the array IWORK. LIWORK >= 1.  If IJOB = 1, 2 or
  233.              4, LIWORK >=  N+6.  If IJOB = 3 or 5, LIWORK >= MAX(2*M*(N-M),
  234.              N+6).
  235.  
  236.              If LIWORK = -1, then a workspace query is assumed; the routine
  237.              only calculates the optimal size of the IWORK array, returns this
  238.              value as the first entry of the IWORK array, and no error message
  239.              related to LIWORK is issued by XERBLA.
  240.  
  241.      INFO    (output) INTEGER
  242.              =0: Successful exit.
  243.              <0: If INFO = -i, the i-th argument had an illegal value.
  244.              =1: Reordering of (A, B) failed because the transformed matrix
  245.              pair (A, B) would be too far from generalized Schur form; the
  246.              problem is very ill-conditioned.  (A, B) may have been partially
  247.              reordered.  If requested, 0 is returned in DIF(*), PL and PR.
  248.  
  249. FFFFUUUURRRRTTTTHHHHEEEERRRR DDDDEEEETTTTAAAAIIIILLLLSSSS
  250.      DTGSEN first collects the selected eigenvalues by computing orthogonal U
  251.      and W that move them to the top left corner of (A, B).  In other words,
  252.      the selected eigenvalues are the eigenvalues of (A11, B11) in:
  253.  
  254.                    U'*(A, B)*W = (A11 A12) (B11 B12) n1
  255.                                  ( 0  A22),( 0  B22) n2
  256.                                    n1  n2    n1  n2
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))
  269.  
  270.  
  271.  
  272.      where N = n1+n2 and U' means the transpose of U. The first n1 columns of
  273.      U and W span the specified pair of left and right eigenspaces (deflating
  274.      subspaces) of (A, B).
  275.  
  276.      If (A, B) has been obtained from the generalized real Schur decomposition
  277.      of a matrix pair (C, D) = Q*(A, B)*Z', then the reordered generalized
  278.      real Schur form of (C, D) is given by
  279.  
  280.               (C, D) = (Q*U)*(U'*(A, B)*W)*(Z*W)',
  281.  
  282.      and the first n1 columns of Q*U and Z*W span the corresponding deflating
  283.      subspaces of (C, D) (Q and Z store Q*U and Z*W, resp.).
  284.  
  285.      Note that if the selected eigenvalue is sufficiently ill-conditioned,
  286.      then its value may differ significantly from its value before reordering.
  287.  
  288.      The reciprocal condition numbers of the left and right eigenspaces
  289.      spanned by the first n1 columns of U and W (or Q*U and Z*W) may be
  290.      returned in DIF(1:2), corresponding to Difu and Difl, resp.
  291.  
  292.      The Difu and Difl are defined as:
  293.  
  294.           Difu[(A11, B11), (A22, B22)] = sigma-min( Zu )
  295.      and
  296.           Difl[(A11, B11), (A22, B22)] = Difu[(A22, B22), (A11, B11)],
  297.  
  298.      where sigma-min(Zu) is the smallest singular value of the (2*n1*n2)-by-
  299.      (2*n1*n2) matrix
  300.  
  301.           Zu = [ kron(In2, A11)  -kron(A22', In1) ]
  302.                [ kron(In2, B11)  -kron(B22', In1) ].
  303.  
  304.      Here, Inx is the identity matrix of size nx and A22' is the transpose of
  305.      A22. kron(X, Y) is the Kronecker product between the matrices X and Y.
  306.  
  307.      When DIF(2) is small, small changes in (A, B) can cause large changes in
  308.      the deflating subspace. An approximate (asymptotic) bound on the maximum
  309.      angular error in the computed deflating subspaces is
  310.  
  311.           EPS * norm((A, B)) / DIF(2),
  312.  
  313.      where EPS is the machine precision.
  314.  
  315.      The reciprocal norm of the projectors on the left and right eigenspaces
  316.      associated with (A11, B11) may be returned in PL and PR.  They are
  317.      computed as follows. First we compute L and R so that P*(A, B)*Q is block
  318.      diagonal, where
  319.  
  320.           P = ( I -L ) n1           Q = ( I R ) n1
  321.               ( 0  I ) n2    and        ( 0 I ) n2
  322.                 n1 n2                    n1 n2
  323.  
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334. DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))
  335.  
  336.  
  337.  
  338.      and (L, R) is the solution to the generalized Sylvester equation
  339.  
  340.           A11*R - L*A22 = -A12
  341.           B11*R - L*B22 = -B12
  342.  
  343.      Then PL = (F-norm(L)**2+1)**(-1/2) and PR = (F-norm(R)**2+1)**(-1/2).  An
  344.      approximate (asymptotic) bound on the average absolute error of the
  345.      selected eigenvalues is
  346.  
  347.           EPS * norm((A, B)) / PL.
  348.  
  349.      There are also global error bounds which valid for perturbations up to a
  350.      certain restriction:  A lower bound (x) on the smallest F-norm(E,F) for
  351.      which an eigenvalue of (A11, B11) may move and coalesce with an
  352.      eigenvalue of (A22, B22) under perturbation (E,F), (i.e. (A + E, B + F),
  353.      is
  354.  
  355.       x = min(Difu,Difl)/((1/(PL*PL)+1/(PR*PR))**(1/2)+2*max(1/PL,1/PR)).
  356.  
  357.      An approximate bound on x can be computed from DIF(1:2), PL and PR.
  358.  
  359.      If y = ( F-norm(E,F) / x) <= 1, the angles between the perturbed (L', R')
  360.      and unperturbed (L, R) left and right deflating subspaces associated with
  361.      the selected cluster in the (1,1)-blocks can be bounded as
  362.  
  363.       max-angle(L, L') <= arctan( y * PL / (1 - y * (1 - PL * PL)**(1/2))
  364.       max-angle(R, R') <= arctan( y * PR / (1 - y * (1 - PR * PR)**(1/2))
  365.  
  366.      See LAPACK User's Guide section 4.11 or the following references for more
  367.      information.
  368.  
  369.      Note that if the default method for computing the Frobenius-norm- based
  370.      estimate DIF is not wanted (see DLATDF), then the parameter IDIFJB (see
  371.      below) should be changed from 3 to 4 (routine DLATDF (IJOB = 2 will be
  372.      used)). See DTGSYL for more details.
  373.  
  374.      Based on contributions by
  375.         Bo Kagstrom and Peter Poromaa, Department of Computing Science,
  376.         Umea University, S-901 87 Umea, Sweden.
  377.  
  378.      References
  379.      ==========
  380.  
  381.      [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
  382.          Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
  383.          M.S. Moonen et al (eds), Linear Algebra for Large Scale and
  384.          Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
  385.  
  386.      [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
  387.          Eigenvalues of a Regular Matrix Pair (A, B) and Condition
  388.          Estimation: Theory, Algorithms and Software,
  389.          Report UMINF - 94.04, Department of Computing Science, Umea
  390.  
  391.  
  392.  
  393.                                                                         PPPPaaaaggggeeee 6666
  394.  
  395.  
  396.  
  397.  
  398.  
  399.  
  400. DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))                                                          DDDDTTTTGGGGSSSSEEEENNNN((((3333SSSS))))
  401.  
  402.  
  403.  
  404.          University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
  405.          Note 87. To appear in Numerical Algorithms, 1996.
  406.  
  407.      [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
  408.          for Solving the Generalized Sylvester Equation and Estimating the
  409.          Separation between Regular Matrix Pairs, Report UMINF - 93.23,
  410.          Department of Computing Science, Umea University, S-901 87 Umea,
  411.          Sweden, December 1993, Revised April 1994, Also as LAPACK Working
  412.          Note 75. To appear in ACM Trans. on Math. Software, Vol 22, No 1,
  413.          1996.
  414.  
  415.  
  416. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  417.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  418.  
  419.      This man page is available only online.
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  
  428.  
  429.  
  430.  
  431.  
  432.  
  433.  
  434.  
  435.  
  436.  
  437.  
  438.  
  439.  
  440.  
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.                                                                         PPPPaaaaggggeeee 7777
  460.  
  461.  
  462.  
  463.